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Abstract 

In this paper we study the performance of a symplectic numerical integrator 
based on the splitting method. This method is applied to a subtle problem i.e. 
higher order resonance of the elastic pendulum. In order to numerically study the 
phase space of the elastic pendulum at higher order resonance, a numerical integrator 
which preserves qualitative features after long integration times is needed. We show 
by means of an example that our symplectic method offers a relatively cheap and 
accurate numerical integrator. 

Keywords. Hamiltonian mechanics, higher-order resonance, elastic pendulum, symplectic nu- 
merical integration, geometric integration. 
AMS classification. 34C15, 37M15, 65P10, 70H08 

1 Introduction 

Higher order resonances are known to have a long time-scale behaviour. From an asymp- 
totic point of view, a first order approximation (such as first order averaging) would not 
be able to clarify the interesting dynamics in such a system. Numerically, this means that 
the integration times needed to capture such behaviour are significantly increased. In this 
paper we present a reasonably cheap method to achieve a qualitatively good result even 
after long integration times. 
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Indonesia 
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Geometric numerical integration methods for (ordinary) differential equations ([0, |T0|, 
131 ) have emerged in the last decade as alternatives to traditional methods (e.g. Runge- 
Kutta methods). Geometric methods are designed to preserve certain properties of a given 
ODE exactly (i.e. without truncation error). The use of geometric methods is particu- 
larly important for long integration times. Examples of geometric integration methods 
include symplectic integrators, volume-preserving integrators, symmetry-preserving inte- 
grators, integrators that preserve first integrals (e.g. energy). Lie-group integrators, etc. 



A survey is given in [IC]. 



It is well known that resonances play an important role in determining the dynamics 
of a given system. In practice, higher order resonances occur more often than lower order 
ones, but their analysis is more complicated. In |12|, Sanders was the first to give an upper 
bound on the size of the resonance domain (the region where interesting dynamics takes 
place) in two degrees of freedom Hamiltonian systems. Numerical studies by van den Broek 
however, provided evidence that the resonance domain is actually much smaller. In 



15| , Tuwankotta and Verhulst derived improved estimates for the size of the resonance 
domain, and provided numerical evidence that for the 4 : 1 and the 6 : 1 resonances of 
the elastic pendulum, their estimates are sharp. The numerical method they used in their 
analysis Q, however, was not powerful enough to be applied to higher order resonances. In 
this paper we construct a symplectic integration method, and use it to show numerically 
that the estimates of the size of the resonance domain in |jl5| are also sharp for the 4 : 3 
and the 3 : 1 resonances. 

Another subtle problem regarding to this resonance manifold is the bifurcation of this 
manifold as the energy increases. To study this problem numerically one would need a 
numerical method which is reasonably cheap and accurate after a long integration times. 

In this paper we will use the elastic pendulum as an example. The elastic pendulum 
is a well known (classical) mechanical problem which has been studied by many authors. 
One of the reasons is that the elastic pendulum can serve as a model for many problems 
in different fields. See the references in |T5|. In itself, the elastic pendulum is a very 
rich dynamical system. For different resonances, it can serve as an example of a chaotic 
system, an auto-parametric excitation system (|T^), or even a linearizable system. The 
system also has (discrete) symmetries which turn out to cause degeneracy in the normal 
form. 

We will first give a brief introduction to the splitting method which is the main ingre- 
dient for the symplectic integrator in this paper. We will then collect the analytical results 
on the elastic pendulum that have been found by various authors. Mostly, in this paper we 
will be concerned with the higher order resonances in the system. All of this will be done 
in the next two sections of the paper. In the fourth section we will compare our symplectic 
integrator with the standard 4-th order Runge-Kutta method and also with an order 7 — 8 
Runge-Kutta method. We end the fourth section by calculating the size of the resonance 
domain of the elastic pendulum at higher order resonance. 



Runge-Kutta method of order 7 — 8 
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2 Symplectic Integration 



Consider a symplectic space = M^" , n G N where each element ^ in has coordinate 
{q,p) and the symplectic form is dq A dp. For any two functions F,G & C°°{Q) define 

which is called the Poisson bracket of F and G. Every function H G C°°(fi) generates a 
(Hamiltonian) vector field defined by {qt, H}, {pi, H}, i = 1, . . . ,n. The dynamics of H is 
then governed by the equations of motion of the form 

Qi ={qi, H] 

Pi ={pi,H}, i = l,... ,n. 

Let X and Y be two Hamiltonian vector fields, defined in Q, associated with Hamiltonians 
Hx and Hy in C°°(f2) respectively. Consider another vector field [-^, ^] which is just the 
commutator of the vector fields X and Y . Then [X, Y] is also a Hamiltonian vector field 
with Hamiltonian H[x,y] = {Hx^Hy}- See for example [0,0, ^ for details. 
We can write the flow of the Hamiltonian vector fields X as 

^x;t = exp(tX) =I + tX + i(tX)2 + + . . . 

(and so does the flow of Y). By the BCH formula, there exists a (formal) Hamiltonian 
vector field Z such that 

z = (X + r) + ^[x, r] + ^ ([X, X, y] + [f, r, x]) + 0{e) (i) 

and exp(tZ) = exp(tX)exp(tF), where [X, X, F] = [X, [X, F]] , and so on. Moreover, 
Yoshida (in [|19]) shows that exp(tX)exp(tF)exp(tX) = exp(tZ), where 

z = (2X + y) + ^ ([F, r, X] - [X, x, r]) + o(t^). (2) 

o 

We note that in terms of the flow, the multiplication of the exponentials above means 
composition of the corresponding flow, i.e. ipY;t ° Vx;t- 

Let r G M be a small positive number and consider a Hamiltonian system with Hamil- 
tonian if(^) = Hx{i) + i^y(^), where ^ G fi, and ^ = X + F. Using (|I]) we have that 
^y;t ° ^x;t is (approximately) the flow of a Hamiltonian system 

^ = (X + r) + ^[x, r] + ^ ([X, X, y] + % r, x]) + 0(r3), 

with Hamiltonian 

H, = Hx + Hy + ^ {Hx, Hy} + ^ {{Hx, Hx, H^} + {H^, H^, Hx}) + 0(r=^). 
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Note that {H, K, F} = {H, {K, F}}. This mean that H - Hr = 0{t) or, in other words 



'PY;T°<fx;T = fx+Y{T)+0{T'^). (3) 

As before and using (|^), we conclude that 

ifx-r o Lfy.^ o Lf^.r = lpx+y{t) + O(r^). (4) 



Suppose that ipx-.r and ipY-.r are numerical integrators of system ^ = X and ^ = Y 
(respectively). We can use symmetric composition (see [H) to improve the accuracy of 
i^x+Y;T- If '^y;t and ipx;T are symplectic, then the composition forms a symplectic numerical 



integrator for X + Y. See [0] for more discussion; also |]IOl for references. If we can split 
H into two (or more) parts which Poisson commute with each other (i.e. the Poisson 
brackets between each pair vanish), then we have H = H^-. This implies that in this case 
the accuracy of the approximation depends only on the accuracy of the integrators for X 
and Y. An example of this case is when we are integrating the Birkhoff normal form of a 
Hamiltonian system. 



3 The Elastic Pendulum 

Consider a spring with spring constant s and length lo to which a mass m is attached. Let 
g be the gravitational constant and / the length of the spring under load in the vertical 
position, and let r be the distance between the mass m and the suspension point. The 
spring can both oscillate in the radial direction and swing like a pendulum. This is called 
the elastic pendulum. See Figure ([I|) for illustration and [|T^ (or [^) for references. 



/ / / / / / / / / / / 



/ / / / / / / / / / / T 



Figure 1: The elastic pendulum. 

The phase space is with canonical coordinate ^ = {z, (p,Pz,Pif), where z = {r — lo)/lo- 
Writing the linear frequencies of the Hamiltonian as Uz = \fsjm and cj^ = \fgjl., the 
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Hamiltonian of the elastic pendulum becomes 




H = ^[P'^ + J7^2) + {-) 1 ---/(- + l)cos^, (5) 

where a = mt^. By choosing the right physical dimensions, we can scale out a. We remark 
that for the elastic pendulum as illustrated in Figure |I], we have uo^ < uj^p. See [|15| for 
details. It is clear that this system possesses symmetry 



T : {z, (p, pz,p^, t) ^ {z, -(f, p,, -p^, t) (6) 
and the reversing symmetries 

Ri : {z,ip,p„p^,t) ^ {z,(p,-p„-p^,-t), 
R2 : {z,^,Pz,p^,t) {z,-ip,-p^,p^,-t). 

If there exist two integers ki and k2 such that kiUz + k2UJ^ = 0, then we say Uz and 
are in resonance. Assuming {\ki\, \ k2\) = 1, we can divide the resonances in two types, e.g. 
lower order resonance if + |A;2| < 5 and higher order resonance if |A;i| + |A;2| > 5. In the 
theory of normal forms, the type of normal form of the Hamiltonian is highly dependent 
on the type of resonance in the system. See |Q. 

In general, the elastic pendulum has at least one fixed point which is the origin of 
phase space. This fixed point is elliptic. For some of the resonances, there is also another 
fixed point which is of the saddle type, i.e. {z, ^p,Pz,P!p) = {—'^iy'ip/^zY, tt, 0, 0). From the 
definition of z, it is clear that the latter fixed point only exists for uozj^tp > V^- The elastic 
pendulum also has a special periodic solution in which (p = p^p = (the normal mode). 
This normal mode is an exact solution of the system derived from (|^). We note that there 
is no nontrivial solution of the form (0, fit), 0,pp(t)). 

Now we turn our attention to the neighborhood of the origin. We refer to for the 
complete derivation of the following Taylor expansion of the Hamiltonian (we have dropped 
the bar) 

H = H2 + eHs + e^H^ + e'^H^ + ■ • • , (8) 

with 

H2 =\ujz {z^ + pI) + \uj^ (<^2 + 

= {rt^'pi - ^^^) 

^5 = -^(^V + 2^.X) 
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In [|T^ the 2 : 1-resonance of the elastic pendulum has been studied intensively. At this 
specific resonance, the system exhibits an interesting phenomenon called auto-parametric 
excitation, e.g. if we start at any initial condition arbitrarily close to the normal mode, 
then we will see energy interchanging between the oscillating and swinging motion. In 
0, the author shows that the normal mode solution (which is the vertical oscillation) is 
unstable and therefore, gives an explanation of the auto-parametric behavior. 

Next we consider two limiting cases of the resonances, i.e. when uj^/oj^p —>■ oo and 
^z/^^p ^ 1- The first limiting case can be interpreted as a case with a very large spring 
constant so that the vertical oscillation can be neglected. The spring pendulum then 
becomes an ordinary pendulum; thus the system is integrable. The other limiting case is 
interpreted as the case where /o = (or very weak spring) 0. Using the transformation r = 
l{z + 1), X = r cos if and y = r sin (f, we transform the Hamiltonian (^ to the Hamiltonian 
of the harmonic oscillator. Thus this case is also integrable. Furthermore, in this case all 
solutions are periodic with the same period which is known as isochronism. This means that 
we can remove the dependence of the period of oscillation of the mathematical pendulum 
on the amplitude, using this specific spring. We note that this isochronism is not derived 
from the normal form (as in [jl8|) but exact. 

All other resonances are higher order resonances. In two degrees of freedom (which is 
the case we consider), for fixed small energy the phase space of the system near the origin 
looks like the phase space of decoupled harmonic oscillator. A consequence of this fact is 
that in the neighborhood of the origin, there is no interaction between the two degrees of 
freedom. The normal mode (if it exists), then becomes elliptic (thus stable). 

Another possible feature of this type of resonance is the existence of a resonance mani- 
fold containing periodic solutions (see paragraph 4.8). We remark that the existence of 
this resonance manifold does not depend on whether the system is integrable or not. In the 
resonance domain (i.e. the neighborhood of the resonant manifold), interesting dynamics 
(in the sense of energy interchanging between the two degrees of freedom) takes place (see 
||12||). Both the size of the domain where the dynamics takes place and the time-scale of 
interaction are characterized by e and the order of the resonance, i.e. the estimate of the 
size of the domain is 

de = [e^) (9) 

and the time-scale of interactions is 0(e^^^) for Uz ■ uj^ = m : n with (m,n) = 1. We 
note that for some of the higher order resonances where uJz/uj^ ^ 1 the resonance manifold 



fails to exist. See 15 for details. 



^This case is unrealistic for the model illustrated in Figure |^. A more realistic mechanical model with 
the same Hamiltonian (|^) can be constructed by only allowing some part of the spring to swing 

''Due to a particular symmetry, some of the lower order resonances become higher order resonances 
([p^). In those cases, (m,n) = 1 need not hold. 
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4 Numerical Studies on the Elastic Pendulum 



One of the aims of this study is to construct a numerical Poincare map {V) for the elastic 
pendulum in higher order resonance. As is explained in the previous section, interest- 
ing dynamics of the higher order resonances takes place in a rather small part of phase 
space. Moreover, the interaction time-scale is also rather long. For these two reasons, we 
need a numerical method which preserves qualitative behavior after a long time of integra- 
tion. Obviously by decreasing the time step of any standard integrator (e.g. Runge-Kutta 
method), we would get a better result. As a consequence however, the actual computation 
time would become prohibitively long. Under these constraints, we would like to propose 
by means of an example that symplectic integrators offer reliable and reasonably cheap 
methods to obtain qualitatively good phase portraits. 

We have selected four of the most prominent higher order resonances in the elastic 
pendulum. For each of the chosen resonances, we derive its corresponding equations of 
motion from (§). This is done because the dependence on the small parameter e is more 
visible there than in (^). Also from the asymptotic analysis point of view, we know that 
(I) truncated to a sufficient degree has enough ingredients of the dynamics of (|]). 

The map V is constructed as follows. We choose the initial values in such a way 
that they all lie in the approximate energy manifold H2 = Eo & M. and in the section 
S = = {z,ip,pz,p,p)\z = 0,pz > 0}. We follow the numerically constructed trajectory 
corresponding to and take the intersection of the trajectory with section S. The in- 
tersection point is defined as V{^o). Starting from 'P(^o) as an initial value, we go on 
integrating and in the same way we find 'P^(^o); and so on. 

The best way of measuring the performance of a numerical integrator is by comparing 
with an exact solution. Due to the presence of the normal mode solution (as an exact 
solution), we can check the performance of the numerical integrator in this way (we will 
do this in section |4.2| ) . Nevertheless, we should remark that none of the nonlinear terms 
play a part in this normal mode solution. Recall that the normal mode is found in the 
invariant manifold {{z, ip,Pz,p^\ip = Pip = 0} and in this manifold the equations of motion 
of (^ are linear. 

Another way of measuring the performance of an integrator is to compare it with 
other methods. One of the best known methods for time integration are the Runge-Kutta 
methods (see P]). We will compare our integrator with a higher order (7-8 order) Runge- 
Kutta method (RK78). The RK78 is based on the method of Runge-Kutta- Felbergh (p!^). 
The advantage of this method is that it provides step-size control. As is indicated by the 
name of the method, to choose the optimal step size it compares the discretizations using 7- 
th order and 8-th order Runge-Kutta methods. A nice discussion on lower order methods of 
this type, can be found in [|1^] pp. 448-454. The coefficients in this method are not uniquely 
determined. For RK78 that we used in this paper, the coefficients were calculated by C. 
Simo from the University of Barcelona. We will also compare the symplectic integrator 
(SI) to the standard 4-th order Runge-Kutta method. 

We will first describe the splitting of the Hamiltonian which is at the core of the 
symplectic integration method in this paper. By combining the fiow of each part of the 
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Hamiltonian, we construct a 4-th order symplectic integrator. The symplecticity is obvious 
since it is the composition of exact Hamiltonian flows. Next we will show the numerical 
comparison between the three integrators, RK78, SI and RK4. We compare them to 
an exact solution. We will also show the performance of the numerical integrators with 
respect to energy preservation. We note that SI are not designed to preserve energy (see 
|T0|). Because RK78 is a higher order method (thus more accurate), we will also compare 
the orbit of RK4 and SI. We will end this section with results on the size of the resonance 
domain calculated by the SI method. 



4.1 The Splitting of the Hamiltonian 

Consider again the expanded Hamiltonian of the elastic pendulum (|]). We split this 
Hamiltonian into integrable parts: H = + + H^, where 

1 2 _ 2I 4 _ ^3^^ 4 

2^ ^ 2V 24^ ^ 



=\uz{z^+pl) + \u^ {v'+pI). 

Note that the equations of motion derived from each part of the Hamiltonian can be 
integrated exactly; thus we know the exact flow v^i.,-, ip2:T) and i^^.^r corresponding to 
H^,H^, and respectively. This splitting has the following advantages. 

• It preserves the Hamiltonian structure of the system. 

• It preserves the symmetry (^) and reversing symmetries (|^) of H. 

• and are of 0{e) compared with H (or H^). 

Note that, for each resonance we will truncate (0) up to and including the degree where 
the resonant terms of the lowest order occur. 
We define 



'^T = <fl;T/2 ° <f2;T/2 « <fS;T O <f2;T/2 « <fl;T/2- 



fir 



From section |^ we know that this is a second order method. Next we define 7 = 1/(2 — v^) 
and ipr = y^-yr ° y^{i-2'y)T ° y^-yr to get a fourth order method. This is known as the generalized 



Yoshida method (see |jTO[)- By, Symplectic Integrator (SI) we will mean this fourth order 
method. This composition preserves the symplectic structure of the system, as well as the 
symmetry (^ and the reversing symmetries (0). This is in contrast with the Runge-Kutta 
methods which only preserves the symmetry (^, but not the symplectic structure, nor the 
reversing symmetries (^). As a consequence the Runge-Kutta methods do not preserve the 
KAM tori caused by symplecticity or reversibility. 
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4.2 Numerical Comparison between RK4, RK78 and SI 



We start by comparing the three numerical methods, i.e. RK4, RK78, and SI. We choose 
the 4 : 1-resonance, which is the most prominent higher order resonance, as a test problem. 
We fix the value of the energy {H2) to be 5 and take e = 0.05. Starting at the initial 
condition z{0) = 0,(p{0) = 0,Pz = a/5/2, and p^{0) = 0, we know that the exact solution 
we are approximating is given by sin(4t),0, v/572 cos(4t),0). We integrate the 

equations of motion up to t = 10^ seconds and keep the result of the last 10 seconds to 
have time series z{tn) and jhitn) produced by each integrator. Then we define a sequence 
Sn = 99990+5r;,/100, n = 0,1, . . . , 200. Using an interpolation method, for each of the time 
series we calculate the numerical In figure ^ we plot the error function — z{sn) 

for each integrator. 




9.999 9.9991 



Figure 2: Plots of the error function z{sn) — z{sn) against time. The upper figure is 
the result of RK4, the middle figure is RK78 and the lower figure is of SI. The time of 
integration is 10^ with a time step for RK4 and SI of 0.025. 

The plots in Figure in |^ clearly indicate the superiority of RK78 compared with the 
other methods (due to the higher order method). The error generated by RK78 is of order 
10~^ for an integration time of 10^ seconds. The minimum time step taken by RK78 is 
0.0228 and the maximum is 0.0238. The error generated by SI on the other hand, is of 
order 10~^. The CPU time of RK78 during this integration is 667.75 seconds. SI completes 
the computation after 446.72 seconds while RK4 only needs 149.83 seconds. 

We will now measure how well these integrators preserve energy. We start integrating 
from an initial condition z{0) = 0, v^(0) = 1.55, Pip{0) = and ^2(0) is determined from 

= 5 (in other word we integrate on the energy manifold H = 5 + 0(e)). The small 
parameter is e = 0.05 and we integrate for t = 10^ seconds. 

For RK78, the integration takes 667.42 second of CPU time. For RK4 and SI we used 
the same time step, that is 10~^. RK4 takes 377.35 seconds while SI takes 807.01 seconds 
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Figure 3: Plots of the energy against time. The sohd hne represents the results from 
SI. The line with '+' represents the results from RK4 and the line with 'o' represents the 
results from RK78. On the left hand plot, we show the results of all three methods with 
the time step 0.01. The time step in the right hand plots is 0.05. The results from RK4 
are plotted seperately since the energy has decreased significantly compared to the other 
two methods. 

of CPU time. It is clear that SI, for this size of time step, is inefficient with regard to CPU 
time. This is due to the fact that to construct a higher order method we have to compose 
the flow several times. We plot the results of the last 10 seconds of the integrations in 
Figure |. We note that in these 10 seconds, the largest time step used by RK78 is 0.02421 
. . . while the smallest is 0.02310 .... It is clear from this, that even though the CPU time 
of RK4 is very good, the result in the sense of conservation of energy is rather poor relative 
to the other methods. 

We increase the time step to 0.05 and integrate the equations of motion starting at the 
same initial condition and for the same time. The CPU time of SI is now 149.74 while for 
the RK4 it is 76.07. Again, in Figure |^ (the right hand plots) we plot the energy against 
time. A significant difference between RK4 and SI then appears in the energy plots. The 
results of symplectic integration are still good compared with the higher order method 
RK78. On the other hand, the results from RK4 are far below the other two. 

4.3 Computation of the Size of the Resonance Domain 

Finally, we calculate the resonance domain for some of the most prominent higher order 
resonances for the elastic pendulum. In Figure ^ we give an example of the resonance 
domain for the 4 : 1 resonance. We note that RK4 fails to produce the section. On the 
other hand, the results from SI are still accurate. We compare the results from SI and 
RK78 in Figure ^. After 5 x 10^ seconds, one loop in the plot is completed. For that time 
of integration, RK78 takes 34.92 seconds of CPU time, while SI takes only 16.35 seconds. 
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Figure 4: Resonance domain for the 4 : 1-resonance. The plots on the left are the results 
from SI while the right hand plots are the results from RK78. The vertical axis is the 
axis and the horizontal axis is (p. The time step is 0.05 and e = 0.05. In the top figures, 
we blow up a part of the pictures underneath. 



This is very useful since to calculate for smaller values of e and higher resonance cases, the 
integration time is a lot longer which makes it impractical to use RK78. 

In Table |l] we list the four most prominent higher order resonances for the elastic 
pendulum. This table is adopted from ||l5l where the authors list six of them. 

The numerical size of the domain in table is computed as follows. We first draw 
several orbits of the Poincare maps V. Using a twist map argument, we can locate the 
resonance domain. By adjusting the initial condition manually, we then approximate the 
heteroclinic cycle of V. See figure ^ for illustration. Using interpolation we construct the 
function ro{6) which represent the distance of a point in the outer cycle to the origin and 6 
is the angle with respect to the positive horizontal axis. We do the same for the inner cycle 
and then calculate max^i \ro{0) — ri{6)\. The higher the resonance is, the more difficult to 
compute the size of the domain in this way. 

For resonances with very high order, manually approximating the heteroclinic cycles 
would become impractical, and one could do the following. First we have to calculate 
the location of the fixed points of the iterated Poincare maps numerically. Then we can 
construct approximations of the stable and unstable manifolds of one of the saddle points. 
By shooting to the next saddle point, we can make corrections to the approximate stable 
and unstable manifold of the fixed point. 
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Resonance 


Resonant 


Analytic 


Numerical 


Error 




part 


loge(4) 


loge(4) 




4 : 1 


H5 


1/2 


0.5091568 


0.004 


6 : 1 


Hj 


3/2 


1.5079998 


0.05 


4 : 3 


H7 


3/2 


1.4478968 


0.09 


3 : 1 


Hs 


2 


2.0898136 


0.35 



Table 1: Comparison between the analytic estimate and the numerical computation of 
the size of the resonance domain of four of the most prominent higher order resonances of 
the elastic pendulum. The second column of this table indicates the part of the expanded 
Hamiltonian in which the lowest order resonant terms are found. 




Figure 5: Plots of log{ds) against log(e) for various resonances. The 4 : 1-resonance is 
plotted using '— o', the 3 : 1-resonance is using ' — h', the 4 : 3-resonance is using '— x' and 
the 6 : 1 resonance is using '— *'. 

5 Discussion 

In this section we summarize the previous sections. First the performance of the integrators 
is summarized in table ^ 

As indicated in table for the 4 : 3 and the 3 : 1 resonances, the higher order 
Runge-Kutta method fails to produce the section. This is caused by the dissipation term, 
artificially introduced by this numerical method, which after a long time of integration 
starts to be more significant. On the other hand, we conclude that the results of our 
symplectic integrator are reliable. This conclusion is also supported by the numerical 
calculations of the size of the resonance domain (listed in Table 

In order to force the higher order Runge-Kutta method to be able to produce the section, 
one could also do the following. Keeping in mind that RK78 has automatic step size control 
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Integrators 


RK4 


RK78 


SI 


The 4 : 1 
resonance 


CPU time {At » 0.025, t = 10^ sec.) 




f\f\/ / ^ QPP 




Preservation oi H (t = lo^ sec.) 


Poor 


Good 


Good 


Orbital Quahty 


Poor 


Very good 


Good 


Section Quahty 




Good 


Good 


The 6 : 1 

resonance 


Orbital Quality 


Poor 


Good 


Good 


Section Quality 




Good 


Good 


The 4 : 3 
resonance 


Orbital Quality 


Poor 


Good 


Good 


Section Quality 






Good 


The 3 : 1 
resonance 


Orbital Quality 


Poor 


Poor 


Good 


Section Quality 






Good 



Table 2: Summary of the performance of the integrators. A bar — indicates that it is not 
feasible to obtain a surface of section for this resonance using this integrator. 



based on the smoothness of the vector field, one could manually set the maximum time step 
for RK78 to be smaller than 0.02310. This would make the integration times extremely 
long however. 

We should remark that in this paper we have made a number of simplifications. One 
is that we have not used the original Hamiltonian. The truncated Taylor expansion of (^) 
is polynomial. Somehow this may have a smoothing effect on the Hamiltonian system. 
It would be interesting to see the effect of this simplification on the dynamics of the full 
system. Another simplification is that, instead of choosing our initial conditions in the 
energy manifold H = C, we are choosing them in = C . By using the full Hamiltonian 
instead of the truncated Taylor expansion of the Hamiltonian, it would become easy to 
choose the initial conditions in the original energy manifold. Nevertheless, since in this 
paper we always start in the section S, we know that we are actually approximating the 
original energy manifold up to order e^. 

We also have not used the presence of the small parameter e in the system. As noted 
in [P], it may be possible to improve our symplectic integrator using this small parameter. 
Still related to this small parameter, one also might ask whether it would be possible to 
go to even smaller values of e. In this paper we took < e < e~°'^. As noted in the 
previous section, the method that we apply in this paper can not be used for computing 
the size of the resonance domain for very high order resonances. This is due to the fact 
that the resonance domain then becomes exceedingly small. This is more or less the same 
difficulty we might encounter if we decrease the value of e. 

Another interesting posibility is to numerically follow the resonance manifold as the 
energy increases. As noted in the introduction, this is numerically difficult problem. Since 
this symplectic integration method offers a cheap and accurate way of producing the res- 
onance domain, it might be posible to numerically study the bifurcation of the resonance 
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manifold as the energy increases. Again, we note that to do so we would have to use the 
full Hamiltonian. 
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